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Abstract. We revisit the effect of peculiar velocities on low-redshift type la supernovae. 
Velocities introduce an additional guaranteed source of correlations between supernova mag¬ 
nitudes that should be considered in all analyses of nearby supernova samples but has largely 
been neglected in the past. Applying a likelihood analysis to the latest compilation of nearby 
supernovae, we find no evidence for the presence of these correlations, although, given the 
significant noise, the data is also consistent with the correlations predicted for the standard 
ACDM model. We then consider the dipolar component of the velocity correlations — the 
frequently studied “bulk velocity” — and explicitly demonstrate that including the velocity 
correlations in the data covariance matrix is crucial for drawing correct and unambiguous 
conclusions about the bulk flow. In particular, current supernova data is consistent with 
no excess bulk flow on top of what is expected in AGDM and effectively captured by the 
covariance. We further clarify the nature of the apparent bulk flow that is inferred when the 
velocity covariance is ignored. We show that a significant fraction of this quantity is expected 
to be noise bias due to uncertainties in supernova magnitudes and not any physical peculiar 
motion. 
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1 Introduction 

Motions of objects in the universe are not entirely random. Objects which are physically close 
to one another respond similarly to the pull of large-scale structure, and as a result their 
peculiar velocities are correlated. Correlations between galaxy peculiar velocities are an old 
subject [1-7], and these velocities have already been used to constrain cosmological models, in 
particular the amount of matter in the universe (see [8] for a review). More recently, peculiar 
velocities have become important in the analysis of type la supernova (SN la) data. At low 
redshift (z < 0.05), typical peculiar velocities of ~300 km/s are a significant contribution to 
the SN redshift (for instance, cz = 3,000 km/s at z = 0.01). These peculiar velocities are a 
nuisance if one is interested in using the SNe to constrain expansion history and dark energy, 
and it is common practice to propagate this extra dispersion into the error budget (e.g. add 
300 km/s X 5/(c2;lnl0) in quadrature to the statistical uncertainty of each SN magnitude). 
However, this neglects significant covariance between the velocities of different SNe. 

Alternatively, one can consider the SN peculiar velocity held itself to be a signal, one 
that should contain useful information about the amount and distribution of matter in the 
universe. Nearby SNe are much fewer in number than nearby galaxies, and given the volume 
limitation for both, this will likely still be the case in the future. On the other hand, SNe are 
more useful on a per-object basis because their individual distances can be inferred directly 
and with relative precision — roughly 7% for each SN, depending on the quality of the 
observations. Therefore, there has been a resurgence of interest in how SN peculiar velocities 
are modeled and used [9, 10]. 

In this paper, we perform a careful study of the SN velocity correlations in current data, 
in particular the way in which they are used to draw conclusions about the so-called “bulk 
velocity” — the motion, relative to the cosmic microwave background (CMB) rest frame, of 
the patch of the universe centered on us and containing the nearby sample. Though we focus 
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on SNe, our methodology is not restricted to SNe and equally applies to analysis of galaxy 
peculiar velocities. 

The paper is organized as follows. In section 2, we review the physics of how peculiar 
velocities affect SN magnitudes. In section 3, we describe the SN samples and how we use 
them. In section 4, we define a general likelihood that is the basis for our analyses, which 
include a test for the presence of velocity correlations (section 5), a test for the presence of 
excess bulk velocity beyond that encoded in the correlations predicted in ACDM (section 6), 
and a comparison to previous work that studied bulk flows without the velocity covariance 
(section 7). We summarize our conclusions in section 8. 


2 Theoretical framework 


2.1 Magnification and SN magnitude residuals at low redshifts 

The magnitude residuals of standard candles like SNe la are directly related to the magnifica¬ 
tion fjL, which is defined as the fractional perturbation in the angular diameter and luminosity 
distances (see [9, 11]), 

1 Adi Ad A ^ ^ 

2 dL{z) dA{z) 

where dA{z) and ^^(z) denote the background distances evaluated at the observed redshift 
z. The second equality, relating luminosity and angular diameter distances, follows from the 
conservation of the photon phase space density. That is, ^ describes both the change in the 
apparent angular size of a spatial ruler as well as the change in observed flux of a standard 
candle. 

Covariant expressions for the magnification at linear order in cosmological perturbations 
have been given in [11-14]. In the conformal-Newtonian (cN) gauge, where the metric is 
written as 

ds"^ = a^(r) [—(1 + 2T) + (1 + 24)) 6ij dx'^dx^] , (2.2) 

the magnification is given, in the notation of [14], by^ 


= 


-2 + 


aHx 


A In a — 2<h 


+ 2k- 2v\\o -If dx{^ -^) 

X Jo 


(2.3) 


where x = x(^) is the coordinate distance inferred using the observed redshift, and 




dxw(x-x)vi (T-4>) 

X 
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Vllo + 



dx [T' - T'] 


(2.4) 

(2.5) 


are the convergence and fractional redshift perturbation, respectively. The latter contains the 
gravitational redshift, Doppler shift, and integrated Sachs-Wolfe effect. Further, denotes 
the Laplacian on a sphere of radius X; denotes the peculiar velocity projected 

along the line of sight n, integrals over x denote integrals along the past lightcone, and a 
subscript o denotes a quantity evaluated at the location of the observer. Note that n, Aina, 

^Here, we have neglected a term that is present if the luminosity of the standard candle depends on time; 
in any case, it is subdominant in the limit we will consider. We have also neglected two pure monopole 
contributions, motivated by the discussion in section 2.2. 
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and any other terms appearing in eq. (2.3) are coordinate-dependent quantities, and only the 
combination given in eq. (2.3) corresponds to an actual (gauge-invariant) observable. This 
can be verified by considering gauge transformations and various test cases [14, 15]. 

For low-redshift SNe, where ^ <C 1 so that x l/(o-f^)) the terms involving the velocity 
are the most significant. This is because the lensing convergence is suppressed for small source 
distances and because the scales probed are much smaller than the horizon. Then, the terms 
involving the potentials and T are also suppressed by roughly aHx relative to the velocity. 
In this case, we obtain 

" -in hi - ’'!») ■ P'6) 

This expression for the magnification, proportional to the relative velocity along the line 
of sight between source and observer, simply arises due to the fact that we evaluate the 
luminosity distance to a supernova using the background distance-redshift relation, while 
the actual redshift is perturbed by the Doppler shift. One can easily verify numerically 
that this approximation is better than 1% for z < 0.1, which is the redshift range we will 
consider in this paper. Note that [9] includes the term —2 in the Aina prefactor in eq. (2.3); 
however, this is not strictly consistent, since the terms involving T in Aina, as well as the 
aberration term —2n||o in eq. (2.3) are of comparable magnitude to this correction (see also 
[16]). We will thus work with eq. (2.6) as the proper low-z limit of eq. (2.3). Note that this 
relation remains valid even if /i becomes of order unity, as long as the velocities ny remain 
small compared to the speed of light. On very small scales [z < 0.01), the velocities are no 
longer described accurately by linear perturbation theory. However, since the SN samples 
considered here are restricted to z > 0.01, we will work with velocities derived from linear 
perturbation theory. Note that, in principle, nonlinear corrections to the velocity could also 
be relevant for higher-redshift SNe, if two SNe happen to be physically close. However, we 
have verihed that nonlinear corrections to the velocities have a small effect (see below). 

As eq. (2.6) shows, the relevant quantity for the magnification at low z is the relative 
velocity between the source and the observer projected along the line of sight. This also 
includes small-scale motions such as the velocity of the Solar System with respect to the 
Milky Way center, which are uncorrelated with large-scale cosmological velocity helds. For 
this reason, we correct the observed SN redshifts to the CMB rest frame using the measured 
CMB dipole moment (see section 3). Then, the magnihcation becomes 

2 

mUcmb = hi - , (2.7) 

where the relevant quantity is now the velocity of the SN relative to the CMB rest frame. 
This simplifies the interpretation, since ny — uy cmb is well described by linear perturbation 
theory. In fact, by performing the calculation in the CMB rest frame (as is normally done), we 
can set vcmb = 0- The following relations will always assume that we work with CMB-frame 
redshifts and in the CMB rest frame. 

It is straightforward to convert a perturbation in the luminosity distance (as in eq. (2.1)) 
into a perturbation of the SN magnitude from the homogeneous background value: 


5m 


5 ^<1 5 V • h 

2 In 10 ^ In 10 aHx 


( 2 . 8 ) 


Note that this relation assumes that /r <C 1 and is thus not applicable at very low redshifts. 
While it is straightforward to derive the proper nonlinear relation for 5m, this is not necessary 
for our purposes, since z > 0.01 in our SN samples. 
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Now consider one object at redshift Zi in direction n* on the sky, and a second at {zj, rij). 
We can derive the covariance of their residuals, 


Sij = {6mi Smj) 


5 

In 10 
5 

In 10 


Q-i Ojj p 

a'iXi a'jXj 

(1 + Zj)"^ (1 + Zj)"^ 

H{zi)dL{zi) H{zj)dL{zj) ’ 


(2.9) 


where ^ij is the velocity covariance given by 

^ij = = ((vi • ni){vj -hj)) 

= ^ ^ ® + '^)ji{kXi)j£{!^Xj)Pe{ni • fij). (2.10) 


Here, primes denote derivatives of the Bessel functions with respect to their arguments, r is 
the conformal time, dr = dt/{a^H), Di is the linear growth function evaluated at redshift 
Zi, and Xi = x{zi). The power spectrum P{k,a) is evaluated in the present (o = 1) and, 
at the large scales we are interested in, only the first ~ 10 terms in the sum over the 
multipoles contribute. As mentioned above, we use velocities derived from linear theory and 
thus insert the linear matter power spectrum for our numerical results. We have verified 
that using a prescription for the nonlinear matter power spectrum in eq. (2.10) does not 
significantly affect our results. We thus conclude that the linear treatment is sufficient for 
our purposes. Physically, this is because the dominant contribution to the covariance comes 
from fairly large-scale modes. Note that in our approach, ((Jm,)^) is assumed to capture the 
random motion contribution to the variance of SN residuals.^ While this is not expected to 
be completely accurate when using the linear matter power spectrum, the difference in the 
diagonal covariance elements is not very significant.^ 

We have denoted this covariance matrix S to emphasize that this is a cosmologically 
guaranteed “signal” to be added to the “noise” covariance matrix that accounts for the 
combination of statistical and systematic errors that affects SN distance measurements, such 
as intrinsic variations in the SN luminosity (see section 3). We again point out that the two 
geometric prefactors in eq. (2.9) each differ by an additive factor of 1 relative to those in 
[9] because we drop the term —2 in eq. (2.3) in order to achieve a consistent low-redshift 
expansion; we have checked that all neglected terms would contribute negligibly at z < 0.1. 


2.2 Monopole subtraction 

The magnification eq. (2.1) and its low-redshift version eq. (2.6) still have a monopole com¬ 
ponent, that is, a contribution that is uniform on the sky. However, since the SN magnitude 
residuals are defined with respect to the best-fit distance-redshift relation, this monopole is 
mostly absorbed in the fit. While there could technically be a residual monopole signal due 

^An alternative approach by [17] models velocities with perturbation theory based on a density field derived 
from other surveys, and complements them with a “thermal” component of ~150 km/s added in quadrature 
to account for nonlinearities. In contrast to our approach, this thus relies on external data sets. A detailed 
comparison between the covariances obtained using these different approaches would be interesting but is 
beyond the scope of this paper. 

^For the low-redshift SNe we consider, nearly all of the redshifts are derived from host galaxy spectra, and 
so the motion of the SN within its host does not contribute to the residuals. 
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to the fact that our fit (to a flat ACDM model, see section 4) is very restricted, we will 
assume here that the bulk of the monopole is removed. Thus, eq. (2.9) needs to be corrected. 
To this end, we define the mean magnitude residual at redshift i as 


Sm{zi) 


6m{zi, n) W (h) , 


( 2 . 11 ) 


where lT(n) is the survey window function, which is normalized such that f Vk(n) d^n = 1. 
Then, noting that we actually measure 6mi = 5mi — 6m{zi), the proper covariance is 

Sij = (dmidniA = ([(5m(zj,nj) — 6m{zi)] [6m{zj,hj) — 5m{zj)\) , (2-12) 


which can be worked out to be 




In 10 


(1 + Zi)"^ (1 + Zj)'^ dDi dDj 


H{zi)dL{zi) H{zj)dL{zj) dr dr 




xEp'’ 


dvr 


Pi{hi ■ Uj) - N(ni) + we{hj)] + dvrlT^ 


(2.13) 


where the survey footprint has been expanded in spherical harmonics, 

lT(n) = '^wimYim{n ), 
im 


and the coefficients Wi{hi) and Wi are defined as 

We{n) = ^^mb"£m(n), Wi = 

m 


(2.14) 


(2.15) 


The extra terms in the square brackets in eq. (2.13) are corrections due to the survey window. 
The Wi are therefore just the angular power spectrum (more precisely, the “pseudo-Q”) of 
the map, while Wi{h) is the i portion of the survey mask at an arbitrary location. Note that, 
due to the required normalization of W, its value where the survey observes is not unity, but 
rather ^ 

I —— (observed sky) 

lT(n) = < ^^sky (2.16) 


(unobserved sky). 


The term in the square parentheses in the last line of eq. (2.13), which includes the 
subtraction of the mean, is therefore a new result that has not, to our knowledge, been 
derived and included in previous analyses (although the existence of such a term has been 
pointed out in [9, 16]). For a full-sky window, it is easy to show that this term becomes 
Pi{ni ■ hj) — 1 for i = 0 and remains equal to the original expression Pi{ni ■ nj) for the other 
multipoles. 

We find that the monopole-subtracted formula leads to small but noticeable changes in 
the results, such as the constraints on the parameter A in section 5, and we recommend that 
it be used in future analyses. 
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3 SN la data and noise covariance 


For our primary SN la dataset, we use the joint light-curve analysis (JLA) [18] of SNe from 
the Supernova Legacy Survey (SNLS) and the Sloan Digital Sky Survey (SDSS). JLA includes 
a recalibration of SNe from the first three years of SNLS [19, 20] along with the complete SN 
sample from SDSS, making it the largest combined SN analysis to date. The final compilation 
includes 740 SNe, ~100 low-redshift SNe from several subsamples, ~350 from SDSS at low 
to intermediate redshifts, ~250 from SNLS at intermediate to high redshifts, and ~10 high- 
redshift SNe observed with the Hubble Space Telescope. 

We combine the individual covariance matrix terms provided^ to compute the full co- 
variance matrix, which includes statistical errors (correlated uncertainties in the light-curve 
measurements, intrinsic scatter, lensing dispersion) and a variety of systematic errors (photo¬ 
metric calibration, uncertainty in the bias correction, light-curve model uncertainty, potential 
non-la contamination, uncertainty in the Milky Way dust extinction correction, and uncer¬ 
tainty in the host galaxy correction). 

Although we compute the covariance matrix as described in the JLA analysis, we 
leave out two contributions to the total error. First, we leave out the additional scatter 
of 150 km/s X 5/(czlnl0) added in quadrature to the other statistical errors on the diago¬ 
nal to account for peculiar velocity. This peculiar velocity scatter does not apply because 
peculiar velocities are not a source of noise in our analysis; instead, they are modeled by the 
formalism discussed in section 2. We also leave out the systematic error term corresponding 
to uncertainty in the peculiar velocity correction applied to the low-z JLA redshifts. Since 
our aim is to study the peculiar velocities themselves, we want to avoid this correction and 
then leave out the systematic error associated with it. To this end, we obtain the CMB-frame 
redshifts zcmb directly from the measured heliocentric redshifts z^ei- Specifically, for each 
SN we compute 


1 + ^CMB = (1 + ^hel) 


1 , t’CMB 

1 H - (ncMB 

c 


n 


(3.1) 


where n is the sky position of the SN, ficMB is the CMB dipole direction, and uqmb is 
the velocity of the Solar System barycenter relative to the CMB rest frame implied by the 
dipole amplitude. We use the measured values ucmb = 369 km/s and iicmb = (1,^) = 
(263.99°,48.26°), where the quoted uncertainties [21] are negligible for our purposes. 

For comparison, we separately consider the Union2 SN la analysis [22] from the Super¬ 
nova Cosmology Project.^ We use the full covariance matrix provided for the Union2 SN 
magnitudes, but as with JLA, we remove the peculiar velocity scatter (300 km/s here) that 
was added to the diagonal. The redshifts given for the Union2 SNe are just the heliocentric 
redshifts transformed to the CMB rest frame, so we use them directly. 

Note that the Union2 compilation of 557 SNe has been superseded by the Union2.1 
compilation [23] of 580 SNe, but here the goal is a fair comparison to previous work that 
analyzes the Union2 data. Since the primary change in Union2.1 is the addition of a set of 
high-redshift SNe, and since only low-redshift SNe are relevant for our analysis, we would 
expect the two compilations to produce very similar results. When substituting Union2.1 for 
Union2, our results do not change qualitatively, but there are some minor differences due to 
new estimates for some corrected SN magnitudes and their errors. Union2.1 also includes a 


^http://supernovae.in2p3.fr/sdss_snls_jla/ 
®http://supernova.Ibl.gov/Union/ 
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host-mass correction (see below) that Union2 does not, but this is relatively small and only 
accounts for part of the magnitude differences. 

Finally, we briefly consider the low-z compilation from the older analysis of [24], also 
for comparison with other work. This compilation does not include an analysis of systematic 
errors, so the uncertainty in a SN magnitude is just a combination of the light-curve measure¬ 
ment errors and the derived intrinsic scatter of (Tint = 0.08 mag. Again, we do not include 
the peculiar velocity scatter of 300 km/s prescribed for a cosmological analysis. Although 
CMB-frame redshifts are given, we transform the given heliocentric redshifts into CMB-frame 
redshifts ourselves using eq. (3.1). 

Because SNe la are not perfect standard candles, it is necessary to correct the observed 
peak magnitude of each SN for the empirical correlations between the SN la absolute mag¬ 
nitude and both the stretch (broadness) and color measure associated with the light-curve 
fitter. More recently, it has become common to fit for a constant offset in the absolute 
magnitude for SNe in high-stellar-mass host galaxies. For JLA, the corrected magnitude is 
therefore given by 

= m + a X (stretch) — /3 x (color) -|- PAM, (3.2) 

where a, (3, and AM are nuisance parameters describing, respectively, the stretch, color, and 
host-mass corrections. The measured P = P{Mj, > IO^^Mq) is the probability that the SN 
occurred in a high-stellar-mass host galaxy. Note that, as mentioned above, Union2 does not 
include this host-mass correction; also, the analysis of [24] uses a different light-curve fitter 
with different (but related) light-curve corrections. 

For each of the three datasets, we fix the SN la nuisance parameters to their best- 
fit values from a fit to the Hubble diagram (for JLA, we perform this fit and correct the 
magnitudes ourselves; for the other datasets, we use precorrected magnitudes). In a proper 
cosmological analysis, one should vary the SN nuisance parameters simultaneously with any 
cosmological parameters. In practice, however, the nuisance parameters are well-constrained 
by the Hubble diagram with little dependence on the cosmological model, so holding them 
fixed should be a good approximation, especially for our purposes here. 

4 Likelihood 

We write the full covariance C as the sum of two contributions, C = S -|- N, where S is the 
signal covariance, dominated by velocities at low z and discussed in section 2, and N is the 
noise covariance, described in section 3. 

Assuming a given cosmological model that allows us to calculate S (eq. (2.9)), the 
optimal way to determine whether the data favors peculiar velocities is to consider evidence 
for the detection of the full signal matrix S. This approach uses more information in the 
data than the search for any particular moment, such as the dipole, of the peculiar velocity 
field. 

Here we would like to detect evidence for coherent departures of supernova magnitudes 
from the mean — that is, clustering. To do this, we introduce a new dimensionless parameter 
A and let S — )• A S, where A = 1 for the hducial model. A = 0 corresponds to the case that 
magnitude residuals are purely due to noise and systematics in the SN la data. We would 
like to test whether A is consistent with one and different from zero. Including the new 
parameter A, the full covariance becomes 

C = AS + N, (4.1) 
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where 0 < ^ < oo. 

On the other hand, allowing for an excess bulk flow component is interesting as well, as 
it can be used to search for signatures beyond the fiducial ACDM model and also allows us to 
compare our results with the existing literature (see section 6). In this case, the magnitude 
residuals are affected by an additional bulk velocity Vbuik (e.g. [9]): 

Am?"''- ^ n,) = - (j^) ' Vb.lt. (4.2) 

where Zi and n* are the redshift and sky position of the SN, while Vbuik is a fixed three- 
dimensional vector. In the quasi-Newtonian picture, Vbuik corresponds to the bulk motion of 
the SN sample; however, in the context of non-standard cosmological models (such as those 
breaking homogeneity or isotropy), this should really be seen as a convenient parametrization 
of the dipole of SN magnitude residuals. 

Putting these ingredients together, we construct a multivariate Gaussian likelihood® 


>C(A, Vbuik) oc 


1 



— Am^C ^Am 
2 


where the elements of the vector Am are 


(4.3) 


(Am)i = mf”’' - - Am®''^'^(vbuik), (4.4) 

where are the observed, corrected magnitudes and AI, Ilm) are the theoretical 

predictions for the background cosmological model (see below). The A4 parameter corre¬ 
sponds to the (unknown) absolute calibration of SNe la; we analytically marginalize over it 
in all analyses (e.g. appendix of [25]). 

We emphasize that, since the covariance depends on the parameter A that we are 
interested in constraining, we need to include a term for the 1 /y^|C| prefactor in addition to 
the usual quantity. Since the covariance is a strictly increasing function of A, neglecting 
the prefactor would lead to the clearly erroneous result that the likelihood is a maximum for 
A —> oo. 

The likelihood in eq. (4.3) is the principal tool we will use for our analyses. In this 
most general form, the likelihood depends on two input quantities (four parameters, since 
the velocity has three components): the normalization A of the signal component of the 
covariance matrix and the excess bulk velocity Vbuik not captured by the velocity covariance. 
Note that, in the fiducial model, A = 1 and Vbuik = 0. 

Throughout our analyses, we assume a flat ACDM model (to = — 1, = 0) with free 

parameters fixed to values consistent with data from Planck [26] and other probes. That is, we 
fix Plm = 0.3, physical matter density = 0.14, physical baryon density = 0.0223, 

scalar spectral index Ug = 0.965, and amplitude of scalar fluctuations Ag = 2.22 x 10“®. The 
corresponding derived value of the Hubble constant is h = 0.683, and that of the amplitude of 
mass fluctuations is cig = 0.79. Within the ACDM model, these parameters are determined 
very precisely using Planck data alone, and we have explicitly checked that modest changes 

®Note that SN flux, or a quantity linearly related to it, might be a better choice for the observable than 
the magnitude, given that we expect the error distribution of the former to be more Gaussian than the latter. 
Nevertheless, this choice should not impact our results, as the fractional errors in flux are not too large, and 
we have explicitly checked that the distribution of the observed magnitudes around the mean is approximately 
Gaussian. Therefore we follow most literature on the subject and work directly with magnitudes. 
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Figure 1. Comparison of the signal (left panel) and noise (right panel) contributions to the full 
covariance matrix for the 111 SNe at z < 0.05 from the JLA compilation. 


in the cosmological model, larger than those allowed by Planck, have a negligible effect on 
our results. We can therefore conclude that adding Planck priors and marginalizing over 
these parameters would not significantly affect our constraints. This is not surprising. The 
background cosmology only affects the monopole of the Supernova magnitudes, and even this 
dependence is weak at very low redshifts (since we marginalize over Ai ). Only a much larger 
change in the parameters would affect the expected pairwise covariance, which we do not 
expect to be able to measure precisely in the first place. 

In figure 1, we compare the noise covariance N to the signal covariance S for our hducial 
cosmology. While the noise contribution is typically larger than the signal, the signal is not 
negligible, and it actually dominates for the lowest-redshift SNe. The noise, unlike the signal, 
becomes effectively smaller as more SNe are used in the analysis, making the signal important 
for the whole redshift range considered (see also the discussion in [27]). 

5 Constraints on the amplitnde of signal covariance 

We first consider whether the data itself shows a preference for the presence of the velocity 
(signal) covariance. Therefore, we explore the constraints on A using the likelihood in eq. (4.3) 
and fixing Vbuik = 0. 

The constraints on the parameter A, which determines the fraction of the velocity 
covariance added to the full covariance C, are shown in figure 2, with the numerical results 
given in table 1. We have adopted a uniform prior on A such that our Bayesian posterior is 
proportional to the likelihood in eq. (4.3). All data choices consistently use the available SNe 
out to z = 0.05. This leaves 111 objects in the JLA analysis and 132 in the Union2 analysis. 
We have explicitly checked that the results are insensitive to the precise redshift cutoff; they 
are driven by the lowest-redshift SNe, and z < 0.05 comfortably captures all of them. 

The solid black curve shows JLA, the most current and rigorously calibrated dataset. 
JLA does not rule out the A = 0 hypothesis; in fact, the likelihood peaks near this value. 
Nevertheless, JLA is fully consistent with the standard value A = 1, with a probability of 
0.07 for A > 1. 
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Figure 2. Constraints on the parameter A that quantifies the amount of velocity correlations {A = 1 
is the standard ACDM value). The JLA data are consistent with A = 1 but do not rule out the 
noise-only hypothesis A = 0. JLA and Union2 give somewhat different constraints, though they are 
not statistically inconsistent. Note that differences remain even after restricting to the rather large 
subset of SNe that they have in common (dashed lines). 


The solid red curve shows the result from the Union2 dataset. While it is noticeably 
different than the JLA result, the two likelihoods are mutually consistent; in particular, A = 1 
is a satisfactory fit to both. Nevertheless, Union2 is different in that it strongly disfavors 
A = 0. 

In order to gain additional insight into the difference between the two datasets, we have 
identihed SNe at z < 0.05 that overlap between the two datasets, a total of 96 objects. 
Performing the analysis on this overlap (dashed lines in hgure 2), we see that the results are 
in better agreement but still somewhat disagree, despite both analyses using the same SN set. 
Part of the reason is that JLA and Union2 determine the magnitudes differently; however, 
even some redshifts do not match. We find a root-mean-square (rms) redshift difference of 
1.4% for the object-to-object comparison of the 96 overlapping SNe, and the largest difference 
is 5%. A further exploration of precisely why the SN redshifts and magnitudes differ is beyond 
the scope of this study, but we have explicitly checked that the difference between JLA and 
Union2 for the overlap set is largely due to differences in the estimated apparent magnitudes, 
not the redshifts.^ 

We also make connection to previous work in [28], where the first cosmological con¬ 
straints from the correlations of SN la peculiar velocities have been obtained. Instead of 
parametrizing the covariance with a multiplicative amplitude, they jointly constrained the 
cosmological parameters Qm and as- Using the earlier SN dataset from [24], they found 
constraints broadly consistent with ACDM values. Using the same dataset, we roughly agree 

^We have checked that the sky positions of JLA and Union2 overlapping SNe do precisely agree. We have 
also found that the subset of these SNe that are also found in the sample of [24] — 40 in total — have redshifts 
that mostly agree very well with JLA, but show larger discrepancies with Union2. 
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with [28], our likelihood favoring j 4 ~ 1.4 with a large uncertainty but effectively ruling out 
A = Q. However, we note that the dataset of [24] includes some SNe with extremely low 
redshifts (as low as z = 0.002), which are not in the Hubble flow and for which the assump¬ 
tion of small residuals in eq. (2.8) breaks down (see also discussion after eq. (2.10)). Objects 
at such extremely low redshifts are mutually separated by distances of a few tens of Mpc; 
their relative velocities are therefore expected to have important nonlinear corrections, in 
addition to the linear relation eq. (2.8) breaking down, and the analysis would have to be 
carefully generalized to take this into account. When we exclude all SNe with z < 0.01 from 
the dataset of [24], the likelihood for A actually looks very similar to the JLA constraints in 
figure 2, favoring H = 0 but still statistically consistent with A = 1. 

To summarize, we find that JLA, the most current and rigorous dataset, does not favor 
the presence of SN velocity covariance guaranteed in the ACDM model, but is nonetheless 
consistent with it. We also find that there is considerable variation in the SN data in terms 
of their constraints on the velocity covariance, and in particular that the optimistic-looking 
results found in [28] were due to some very-low-redshift SNe that may be too nearby for 
accurate modeling with linear theory. 

The covariance of SN flux residuals has the potential to provide additional information 
about cosmological parameters and other interesting physics [29-33]. Our results suggest 
that current data do not yield interesting constraints. This will likely change with larger, 
homogeneous samples with greater sky coverage, such as those expected from the Large 
Synoptic Survey Telescope (LSST), currently under construction. 

6 Constraints on excess bnlk velocity 

The optimal way to search for the effect of SN velocity correlations in the context of a fidu¬ 
cial cosmological model is to test for the presence of the full signal covariance, as we did in 
section 5. However, we can also use the SN magnitude residuals to search for a dipolar dis¬ 
tortion (for example, due to bulk motion) beyond what is expected in ACDM. This provides 
constraints on physics beyond the concordance cosmology, such as a breaking of statistical 
isotropy or homogeneity, or the presence of a single large long-wavelength perturbation. The 
search for bulk flows is the subject of a significant body of literature. Here we use a different 
approach, and we make the connection to previous literature in the next section. 

Bulk velocity is usually defined as the motion of the volume spanned by SNe la and the 
rest frame defined by the CMB. Moving the SN redshifts to the CMB frame, we are looking 
for an overall motion between the SN volume and the rest frame. Since we aim to search for 
an excess bulk flow beyond ACDM, we include the velocity correlations S in the likelihood, 
fixing A = 1. Since S includes all velocity correlations within ACDM, including the dipole, 
we expect the posterior for the bulk flow to be consistent with zero if ACDM provides a 
good description of the SN data. Of course, this assumes that our linear modeling of the 
velocity correlations is accurate for the SN sample and also that there are no unaccounted-for 
systematic errors in the SN data that could masquerade as an excess bulk flow. 

For this analysis, we therefore adopt the likelihood from eq. (4.3), setting A = 1 but 
allowing the bulk velocity Vbuik = ('^^buik) 4>) to vary in magnitude and direction. In order to 

get a sense of the three-dimensional likelihood C{vhuik, G, (j^), we first consider the likelihood 
of the excess bulk flow amplitude in a cut through the best-fit direction, that is, as a function 
of Ubuik with 9 and (p set to their maximum-likelihood values. This is shown in figure 3, and 
note that we continue the scan past zero velocity in the direction opposite that of the best 
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Figure 3. A slice through the 3D likelihood for excess bulk velocity. The direction is fixed to be 
fimax-iike (different for each dataset), while the amplitude of the dipole is varied and allowed to be 
positive or negative. Conclusions about the bulk flow would differ significantly if the velocity signal 
covariance were set to zero (dashed lines), as in most previous work on the subject. 

fit by letting the amplitude of the bulk flow take negative values. Because the likelihood 
is Gaussian, and because Vbuik enters the observable magnitude linearly (see eq. (4.2)), the 
likelihood of the bulk velocity is also guaranteed to be Gaussian. Therefore, the likelihood 
ratio between the best-fit (upuik, 0, (j)) and Ubuik = 0 immediately translates into the confidence 
at which zero excess bulk velocity is ruled out, assuming uniform priors on each component 
of the vector Vbuik- 

Figure 3 shows that, once the velocity covariance is properly taken into account, SN 
data do not favor any bulk velocity beyond the amount expected in ACDM. For example, 
the JLA likelihood peaks at (ubuikj^,^) = (187 km/s, 323°, 25°), but this likelihood is larger 
than that for Ubuik = 0 only by —2Aln£ = —1.6, far too low even for 68.3% (Icr) evidence, 
which in three dimensions would be —2Aln/i ~ —3.5. 

We would now like to explicitly calculate the posterior distribution of the amplitude of 
an excess bulk velocity. Assuming uniform priors on each component of Vbuik would produce 
an additional factor in the posterior, driving it to zero for the Ubuik = 0 case so that 
Vbuik = 0 is automatically ruled out. An alternative that allows us to test the Vbuik = 0 
assumption, implicitly (or explicitly [34]) adopted by most previous work, is to choose the 
prior Pr(ubuik) oc 1/u^uib or, alternatively, to consider the angle-averaged likelihood 

pangle-averaged(^^^j^) (X J £(vbulk) ^ COS 0 # . (6.1) 

We plot this likelihood in figure 4. It is immediately apparent that both JLA and 
Union2 data show no preference for excess bulk velocity, though there is a large uncertainty. 
Performing the same analysis but setting the velocity correlations to zero (so that C = N 
with S = 0), the results are drastically different, favoring bulk flows of several hundred km/s 
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Figure 4. Angle-averaged posterior on the amplitude of bulk velocity in excess of the correlations 
captured by the full ACDM velocity covariance. Both JLA and Union2 data are consistent with 
zero velocity, with relatively large error. The conclusion would again be very different if the velocity 
covariance were artificially set to zero (dashed lines). 


and, in the case of Union2, firmly ruling out the Ubuik = 0 case. This is in agreement with the 
conclusion found in previous work [35-40]. In the dashed curves shown in figure 4, we have 
not inclnded the 150 km/s (300 km/s) scatter that is added in quadrature to the diagonal 
of the noise covariance N for JLA (Union2) data in some analyses. This clearly does not 
capture the significant fnll covariance of SN residuals due to large-scale structure. However, 
when adding back this contribntion, we find nearly perfect agreement with the resnlts of [37] . 

Given the importance of this issue, we stress again that the dashed lines in figures 3 
and 4 do not show the proper likelihood of any peculiar velocity, ACDM or otherwise, since 
a gnaranteed component of the covariance of the data has been neglected (nevertheless, one 
can still use these incorrect likelihoods to construct an estimator for bulk flow, which we will 
turn to in the next section). Since the velocity covariance S gives a guaranteed sonrce of 
correlations in ACDM, we argue that it should always be included in likelihood analyses of 
SN magnitude residuals. Neglecting this covariance will lead to suboptimal estimators and, 
in general, biased results. 

7 Relation to previous bulk flow measurements 

A significant body of earlier work on SN velocities neglected the velocity covariance (and 
the lensing covariance, which is important at higher redshifts). These analyses often fonnd 
evidence for non-zero bulk flow, and we confirm these findings with our dashed curves in 
figures 3 and 4. This bulk velocity with A = 0 is difficult to interpret, since it was obtained 
in an analysis with a guaranteed contribution to the covariance of the observables artificially 
set to zero. 
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Nevertheless, one can derive a theoretical expectation for what one should expect for 
the bulk velocity derived in this way; we call it and in some previous work it 

corresponds to what the authors simply call “bulk velocity.” Note that this is not the excess 
bulk velocity over the ACDM expectation considered in the previous sections. 

We will adopt the likelihood from eq. (4.3) once more, but we define a new vector x via 
= ^^buik X, where 


Xi 


5 ^ (1 + Zj)'^ ^ ^ 

Inioy H{zi)dL{zi) 


(7.1) 


To estimate we set the signal covariance to zero, assume a fixed direction 

Abuik, and find the maximum of the likelihood; that is 


'-'noise-only , 

Wulk ^ 


max [C{A = 0, Vbuik)] • 

■''bulk 


(7.2) 


Maximizing with respect to Ubuikj one finds 


'-'noise-only 

^bulk 


yT (5m 
yTx 


(7.3) 


where we have defined for convenience a new vector y as y = N'^^x, and where (5m is the 
vector of SN magnitude residuals (eq. (4.2) without the term). Assuming ACDM, 

and still keeping the direction fibuik fixed, the expected (mean) velocity is of course zero since 

{5mi) = 0. 

We are more interested in the variance of this quantity, which can be computed directly 
and is equal to 


/ 7'.noise-only\2\ _ yT(S + N) y 

y(Wulk ) j (yTx) 2 


(7.4) 


since the true SN magnitude covariance for ACDM is the sum of both signal and noise: 
{{6mi){6mj)) = Sij + Nij. The square root of this quantity is the desired theoretical expec¬ 
tation for the rms bulk velocity in ACDM when one ignores the signal covariance matrix. 

Using the JLA SNe up to redshift z = 0.05, the ACDM expectation for the rms velocity 
varies from about 150-170 km/s as a function of fibuiki with a sky-averaged value of 162 km/s. 
Assuming only noise in the data, {{6mi){5mj)) = Nij, the result is 71 km/s. The Union2 
data give similar results. 

We have therefore found that the predicted rms value of assuming ACDM 

and SN data up to z = 0.05, is ~160 km/s, and that nearly half of this value would be 
generated by statistical scatter in SN magnitudes in the absence of any peculiar velocities in 
the universe (such a contribution is sometimes referred to as the “noise bias”). The peak of 
the JLA likelihood (black dashed line in figure 4) is in agreement with the ACDM expectation; 
Union2 gives a somewhat larger result. Again, however, we caution that the analysis in this 
section is suboptimal, given that we do condense the data into a weighted dipole estimator 
Wuik^"°”*^ rather than using the full covariance. Moreover, this estimator is significantly 
affected by noise bias which needs to be subtracted. Given the uncertainties in the noise 
covariance, the subtraction of noise bias will lead to additional systematic uncertainties in 
the actual peculiar velocity measurement. 

We have not attempted to repeat the analyses of some past work that studied the 
velocity field of low-redshift SNe [38-45] or the anisotropy of the universe from the distribution 
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^bulk 

(fig 3) 

”* (fig 4) 

Survey 

ML 

95% C.L. 

o 

II 

(M 

<1 


ML 

95% C.L. 

95% C.L. 

JLA 

0.19 

(0,1.15) 

0.24 

4.13 

187 

(-108,485) 

(0,376) 

Union2 

1.19 

(0.19,3.27) 

13.2 

0.07 

265 

(-37, 568) 

(0,456) 


Table 1. Summary of numerical results. For both JLA and Union2, we show the best-fit, maximum- 
likelihood (ML) values and 95% confidence intervals for A. We also show the quantity (that is, 
—2AlnL) between the best-fit value and special values A = 0 (no velocity signal) and A = 1 (ACDM 
velocity signal). We also show ML and 95% intervals for bulk velocity in the best-fit direction and 
95% intervals for angle-averaged bulk velocity (here we do not report ML values, which are near zero). 
All velocities are in units of km/s. 


of nearby SN distances [46-51], since these studies adopted a wide variety of approaches and, 
in some cases, complicated statistical procedures whose results are calibrated on simulations. 
We emphasize, however, that the velocity covariance due to large-scale structure should be 
included in any such analyses in order to obtain unbiased results and draw reliable statistical 
conclusions about the velocity held of the nearby universe. 

8 Conclusions 

In this paper we have revisited the constraints on bulk velocity — the relative motion of the 
volume populated by nearby SNe la and the rest frame dehned by the CMB. Our emphasis 
was on a precise and clear procedure for selecting the data, performing the analysis, and 
modeling the theoretical expectation. We concentrated on SNe la as tracers of cosmic struc¬ 
ture and studied two separate (but overlapping) datasets. Our methodology applies equally 
well to galaxies and other tracers of large-scale structure. 

We argued, and demonstrated with explicit calculations, that inclusion of the “signal” 
covariance matrix that captures the peculiar velocity correlations between SNe is crucial. 
The velocities provide a guaranteed source of covariance between SNe; while the velocity 
contribution is subdominant compared to the noise except at the lowest redshifts (see hg- 
ure 1), it does not become smaller as more SNe are included in the analysis. Neglecting the 
velocity covariance, as done by a significant body of earlier work on SN velocities and tests 
of statistical isotropy, leads to results that are both biased and difficult to interpret. 

Our approach was based on a likelihood that includes both the signal and noise co- 
variance and four free parameters: the normalization A, specifying the fraction of the signal 
added to the covariance, and three components of an excess bulk velocity Vbuik beyond that 
which is encoded in the signal covariance. For the fiducial ACDM model, A = 1 and Vbuik = 0. 

We hrst investigated whether the standard ACDM velocity covariance (A = 1) is pre¬ 
ferred over the case in which the covariance is ignored {A = 0); that is, we constrained 
£(A, Vbuik = 0). We found that the JLA dataset, while consistent with A = 1, cannot rule 
out A = 0, and in fact its likelihood peaks near zero (figure 2). Therefore, we did not find 
convincing evidence in the data for the correlations expected from velocities. Although we 
expect things to change with future data, when precise measurements of a quantity like A 
will effectively constrain cosmological parameters such as Dm and ag, our results indicate 
that current data are not close to providing such useful constraints. 
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We then pursued a different approach; we assumed a standard ACDM velocity covari¬ 
ance (the A = 1 case) and tested for excess bulk velocity Vbuik beyond that already captured 
by the covariance — that is, we adopted the likelihood £(A = 1, Vbuik) in the analysis. We 
found that current SN data provide no evidence for a departure from the null hypothesis 
Vbuik = 0 (figures 3 and 4). This result is in sharp contrast to the inference one would have 
made by ignoring the velocity covariance (that is, setting A = 0 in the same analysis), as 
some previous analyses in the literature have done. 

To better understand this latter case, we separately studied inferred constraints on a 
“non-excess” bulk velocity where the velocity covariance has been ignored — that is, the 
C{A = 0, case. Note that this bulk velocity is more difficult to interpret since it 

was obtained by ignoring a guaranteed source of correlations in the data. We showed that the 
rms of the estimated assuming ACDM and SN data up to 2 ; = 0.05, is expected to 

be ~160 km/s, and that nearly half of this value would be generated by a contribution purely 
from intrinsic and observational scatter in the SN magnitudes. Therefore, there are really 
two problems with this approach: not only is the constrained quantity difficult to interpret, 
but it is also guaranteed to be nonzero even without any peculiar velocities in the universe, 
which is clearly not optimal for cosmological interpretations. 

The mapping of velocities in the universe using nearby tracers of large-scale structure 
has had a remarkably long and productive history. With upcoming large-held, fast-scanning 
surveys, it is likely that data will become of sufficiently high quality to enable peculiar veloc¬ 
ities to progress to the next level and become competitive cosmological probes. Of course, 
data analysis and theoretical modeling will have to progress as well. In this paper we have 
demonstrated that, even for current data, clearly dehning the quantities to be constrained 
and carefully accounting for the guaranteed correlations between objects due to large-scale 
structure are two factors of key importance. 
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